.

Séries Temporais

Igor Kuivjogi / Kaíque Ferreira / Nilo Célio

4° Trabalho de Séries Temporais

22 de Agosto de 2020

.

Referências:

  • Time Series Analysis (with Applications in R) - Jonathan D. Cryer, Kung-Sik Chan

Lendo dados

Gráfico da série

ggplotly(
ggplot(renner, aes(x = ref.date, y = ret.adjusted.prices)) +
  geom_line() +
  labs(x = 'data', y = 'retorno (em %)')
)

ACF

# guardando série de retorno
rt <- renner$ret.adjusted.prices

acf(rt)

PACF

pacf(rt)

Os gráficos indicam que a śerie \(r_t\) não possui autocorrelação.
Porém, pelo gráfico da série, podemos ver que a variância não parece ser constante ao longo do tempo, pois em alguns dias seguidos o retorno é muito alto ou muito baixo.
As séries correlacionadas podem ser estudadas usando o valor absoluto ou o valor quadrático. Se os retornos diários são i.i.d., então os retornos diários em valor absoluto (ou ao quadrado) também serão.

ACF para valor absoluto e quadrático

par(mfrow = c(2, 1))
acf(abs(rt), main = 'retorno em valor absoluto')
acf(rt ^ 2, main = 'retorno quadrático')

par(mfrow = c(1, 1))

PACF para valor absoluto e quadrático

par(mfrow = c(2, 1))
pacf(abs(rt), main = 'retorno em valor absoluto')
pacf(rt ^ 2, main = 'retorno quadrático')

par(mfrow = c(1, 1))

Como podemos ver, as séries em valor absoluto e ao quadrado possuem autocorrelação, contradizendo a hipótese (para a série original \(r_t\)) que não existia autocorrelação.

Modelo ARCH

Os modelos ARCH, ou modelos auto-regressivos com heterocedasticidade condicional, foram introduzidos por Engle (1982), com o objetivo de estimar a variância da inflação. A expectativa aqui é que o retorno \(X_t\) seja não-correlacionado serialmente, mas a volatilidade (variância condicional) depende de retornos passados por meio de uma função quadrática.

Definição

Um modelo \(ARCH(p)\) é definido por

\(X_t = \sqrt{h_t}\epsilon_t\),

\(h_t = \alpha_0 + \alpha_1X_{t-1}^2+\ldots+ \alpha_rX_{t-p}^2\),

onde \(\epsilon_t\) é uma sequência de variáveis aleatórias independentes e identicamente distribuídas (i.i.d.) com média zero e variância um, \(\alpha_0 > 0\), \(\alpha_1 \geq 0\), \(i > 0\).

Na prática, usualmente supomos \(\epsilon_t\) ~ \(N(0,1)\) ou \(\epsilon_t\) ~ \(t_v\) (distribuição \(t\) de Student com \(v\) graus de liberdade).

Propriedades do modelo \(ARCH\)

\[ E(X_t) = 0 \]

\[ Var(X_t) = \frac{\alpha_0}{1-\alpha_1} \]

\[ Cov(X_{t+k}X_t|F_{t+k-1})=0 \]

O modelo ARCH na verdade é um modelo de regressão, em que a variável resposta é a volatilidade condicional e os lags passados de \(r_t ^ 2\) são as covariáveis.

O modelo \(ARCH(1)\) considera somente o retorno quadrático no instante \(t\) como função do retorno quadrático anterior, ou seja, \(r_{t-1} ^ 2\). Portanto num modelo \(ARCH(q)\), considera-se os \(t - q\) retornos quadráticos, e não somente o anterior.

Modelo GARCH

Uma generalização dos modelos \(ARCH\) foi sugerida por Bollerslev (1986), chamado modelo \(GARCH\) (generalized \(ARCH\)).

Definição

Um modelo \(GARCH(p,q)\) é definido por

\(X_t = \sqrt{h_t}\epsilon_t\),

\(h_t = \alpha_0 + \sum_{\substack {i=1}}^{r} \alpha_iX_{\substack {t-i}}^{2} + \sum_{\substack {j=1}}^{s} \beta_jX_{\substack {t-j}}\) ,

em que \(\epsilon_t\) i.i.d., \(\alpha_0 > 0\), \(\alpha_i \geq 0\), \(\beta_j \geq 0\), \(\sum_{\substack {i=1}}^{q}(\alpha_i + \beta_i) <1\), \(r = max(p,q)\).

No modelo \(GARCH(p, q)\), além de considerarmos os retornos quadráticos passados, também leva-se em consideração as variâncias condicionais passadas.

Modelo ARMA para série \(r_t\)

O modelo \(GARCH(p, q)\) para a série original de retornos \(r_t\) implica que a série quadrática \(r_t ^ 2\) é modelada por um \(ARMA(max\{p, q\}, q)\). Então podemos usar todas técnicas de identificação que aprendemos até agora para modelar \(r_t^2\) como um modelo \(ARMA(p, q)\):

# deixando de fora 10 valores
rt_treino <- rt[1:(length(rt) - 10)]

# ajustando modelo ARMA
mod <- auto.arima(rt_treino ^ 2, d = 0)
summary(mod)
## Series: rt_treino^2 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ma1    mean
##       0.9838  -0.9677  3.9145
## s.e.  0.0127   0.0174  0.3853
## 
## sigma^2 estimated as 47.44:  log likelihood=-4127.46
## AIC=8262.93   AICc=8262.96   BIC=8283.4
## 
## Training set error measures:
##                       ME     RMSE      MAE  MPE MAPE      MASE       ACF1
## Training set 0.009930523 6.879247 4.040285 -Inf  Inf 0.7831633 0.01974922

O método automático auto.arima sugere um \(ARMA(1, 1)\) para ajustar o retorno quadrático \(r_t ^ 2\).

Diagnóstico do modelo ARMA, para diferentes lags

checkresiduals(mod, plot = F)
## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 9.3029, df = 7, p-value = 0.2316
## 
## Model df: 3.   Total lags used: 10
checkresiduals(mod, plot = F, lag = 12)
## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 9.8054, df = 9, p-value = 0.3665
## 
## Model df: 3.   Total lags used: 12
checkresiduals(mod, plot = F, lag = 16)
## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 17.202, df = 13, p-value = 0.1902
## 
## Model df: 3.   Total lags used: 16
checkresiduals(mod, plot = F, lag = 20)
## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 21.009, df = 17, p-value = 0.2259
## 
## Model df: 3.   Total lags used: 20
checkresiduals(mod, plot = F, lag = 24)
## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 23.412, df = 21, p-value = 0.3224
## 
## Model df: 3.   Total lags used: 24

Para diferentes lags, o teste de Ljung-Box não rejeita a hipótese nula de que os resíduos do modelo \(ARMA(1, 1)\) são ruído branco. Através desta identificação, podemos ajustar um \(GARCH(1, 1)\).

Ajustando o modelo GARCH

g_model <- garch(rt_treino, order = c(1, 1), trace = F)
summary(g_model)
## 
## Call:
## garch(x = rt_treino, order = c(1, 1), trace = F)
## 
## Model:
## GARCH(1,1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.14480 -0.55946  0.02364  0.69758  5.45180 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)    
## a0   0.08554     0.07141    1.198   0.2310    
## a1   0.01682     0.00765    2.199   0.0279 *  
## b1   0.96151     0.02432   39.541   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 78.439, df = 2, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 0.047388, df = 1, p-value = 0.8277

Diagnósticos do modelo

O livro Time Series Analysis (with Applications in R) - Jonathan D. Cryer, Kung-Sik Chan sugere que se o modelo \(GARCH(p, q)\), foi escolhido adequadamente, então o valor absoluto, ou quadrado dos resíduos, padronizados do modelo não devem possuir pontos significativos para o gráfico de ACF:

par(mfrow = c(2, 1))
acf(abs(residuals(g_model)), na.action = na.pass)
acf(residuals(g_model) ^ 2, na.action = na.pass)

par(mfrow = c(1, 1))

Os gráficos sugerem que o modelo \(GARCH(1, 1)\) é adequado aos dados de retorno de ações da Renner.

Função matemática para o modelo \(GARCH(1, 1)\):

\(X_t = \sqrt{h_t}\epsilon_t\),

\(h_t = 0.08554 + \sum_{\substack {i=1}}^{q} 0.01682X_{\substack {t-i}}^{2} + \sum_{\substack {j=1}}^{s} 0.96151X_{\substack {t-j}}\) ,

em que \(\epsilon_t\) i.i.d., \(q = max(p,q)\).

Ajustamento do modelo \(GARCH(1, 1)\) à série \(r_t\)

# eleva ao quadrado porque o modelo ajusta o desvio padrão condicional, e não a variância condicional
ggplotly(
tibble(t = renner$ref.date[1:1233],
       observado = rt_treino,
       GARCH_pos = fitted(g_model)[, 1] ^ 2,
       GARCH_neg = -fitted(g_model)[, 2] ^ 2) %>% 
  ggplot(aes(x = t)) +
  geom_line(aes(y = observado), color = 'grey50') +
  geom_line(aes(y = GARCH_pos), color = 'red') +
  geom_line(aes(y = GARCH_neg), color = 'red4')
)

Gráfico da predição

A predição retorna o desvio padrão condicional, por isso eleva-se ao quadrado.

# dados de teste
teste <- rt[1234:1243]
pred <- predict(g_model, newdata = teste)[, 1] ^ 2

tab_pred <- tibble(
  t = renner$ref.date[1:1243],
  observado = rt[1:1243],
  predito_pos = c(rep(NA, 1233), pred),
  predito_neg = c(rep(NA, 1233), -pred)
)

ggplotly(
tab_pred[1100:1243, ] %>%
  ggplot(aes(x = t)) +
  geom_line(aes(y = observado), color = 'grey50') +
  geom_line(aes(y = predito_pos), color = 'red') +
  geom_line(aes(y = predito_neg), color = 'red4') +
  geom_hline(yintercept = 0, linetype = 'dashed')
)

Erro de predição

rmse <- function(y, y_hat) {
  sqrt(sum((y - y_hat) ^ 2) / length(y))
}

cat('RMSE para dados de treino:\n')
## RMSE para dados de treino:
rmse(rt[1235:1243], predict(g_model, newdata = rt[1234:1243])[, 1][2:10])
## [1] 1.773664
cat('RMSE para dados de teste:\n')
## RMSE para dados de teste:
rmse(teste[2:10], pred[2:10])
## [1] 3.424494